Phase equilibria in the polydisperse Zwanzig model of hard rods 
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We study the phase behaviour of the Zwanzig model of suspensions of hard rods, allowing for 
polydispersity in the lengths of the rods. In spite of the simplified nature of the model (rods are 
restricted to lie along one of three orthogonal axes) , the results agree qualitatively with experimental 
observations: the coexistence region broadens significantly as the polydispersity increases, and strong 
fractionation occurs, with long rods found preferentially in the nematic phase. These conclusions 
are obtained from an analysis of the exact phase equilibrium equations. In the second part of the 
paper, we consider the application of the recently developed "moment free energy method" to the 
polydisperse Zwanzig model. Even though the model contains non-conserved densities due to the 
orientational degrees of freedom, most of the exactness statements (regarding the onset of phase 
coexistence, spinodals, and critical points) derived previously for systems with conserved densities 
remain valid. The accuracy of the results from the moment free energy increases as more and more 
additional moments are retained in the description. We show how this increase in accuracy can be 
monitored without relying on knowledge of the exact results, and discuss an adaptive technique for 
choosing the extra moments optimally. 



I. INTRODUCTION 

Solutions of rod-like particles undergo a transition 
from a disordered isotropic phase to an ordered nematic 
phase as the concentration of rods is increased. In the 
isotropic phase the rods have no preferred orientation, 
whereas in the nematic phase there is a favoured average 
alignment of the_rods. Such a transition has been ob- 
served by Zocheiu in solutions of rod-like V2O5 particles, 
and later by Bawden, Pirie, Bernal and FankuchennB in. 
solutions of tobacco mosaic virus. In 1949, Onsagero 
showed that the ordering could be explained, at least 
for solutions of monodisperse long thin rigid rods, by 
considering only the competition between the excluded 
volume interaction and the orientational entropy. The 
key to his model is the distribution function, P{6), that 
represents the fraction of rods with a given angular ori- 
entation, 9, with respect to some director. In principle, 
this function can be determined self-consistently by min- 
imizing the free energy; however, the resultant non-linear 
integral equation for P{9) cannot be solved analytically. 
To progress further, Onsager introduced a trial function 
with a single variational parameter, a. By now mini- 
mizing the resultant free energy with respect to a, its 
concentration dependence can be calculated. From this 
dependence the conditions for coexistence of the isotropic 
and nematic phases can be determined. 

There are several difficulties in quantitatively com- 
paring experimental results, from systems such as solu- 
tions of tobacco mosaic virus, with Onsager 's predictions. 
Firstly, in addition to the hard core repulsion, the inter- 
action between rods often includes a soft repulsion due to 



electrostatic forces. Secondly, the rods are rarely rigid, 
and semi-flexibility must be accounted for. To overcome 
these difficulties, Buining, Veldhiiizen, Pathmamanoha- 
ran, Hansen and Lekkerkerkeidij synthesized stericallK 
stabilized beohmite particles. Buining and LekkerkerkeiCI 
and van Bruggen, van der Kooij and Lekkerkerkem then 
studied the phase behaviour of solutions of these parti- 
cles, which can be modelled reasonably with a rigid hard 
rod interaction, with corrections arising from soft elec- 
trostatic repulsions and semi-flexibility being negligible. 
However, there remains one further complication: as with 
most polymeric systems, the particles are polydisperse. 

The Onsager approach has_ be/an modified to |eiiable 
phase diagrams for bidispersm and tridispersetd sys- 
tems, in which the rods differ only in length, to be cal- 
culated. A rich variety of behaviour has been predicted, 
such as widening of the region of coexistence, fractiona- 
tion of the longer rods into the nematic phase, nematic- 
nematic coexistence, and re-entrant phases. Qualita- 
tively, the first three of these are in agreement with the 
experiments reported in Ref. ^ and ^, and the widening 
of the region of coexistence between phases appears to 
be a general feature of polydisperse systems. Despite the 
obvious difficulties in trying to map a continuous poly- 
disjjerse system onto a bidisperse mixture, Merchant and 
Rillt3 have attempted to analyse the transition concen- 
tration in solutions of polydisperse rod- like DNA in terms 
of the theory presented in Ref. It is clear that quanti- 
tative agreement is still lacking. The generalisation of the 
Onsager approach to continuous polydispersity is, how- 
ever, far from simple. The only attempts, to date, that 
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we are aware of, have treated the polydispersity pertur- 
bativelytSEj; this hmits the vahdity of the analysis to sit- 
uations with rather narrow distributions of rod lengths. 

Another approach to the problem of niopodisperse rod- 
like mixtures was adopted by Zwanzig)l3: he restricted 
the orientations of the rods to be in one of three mu- 
tually perpendicular directions. This enables the exact 
calculation of higher order virial coefficients, and the ori- 
entational distribution can be determined without ap- 
proximations. In contrast to the Onsager approach, the 
Zwanzig model may be readily extended to polydisperse 
systems. For bidisperse mixtures it has already been 
shown by Clarke and McLeishE] that the qualitative fea- 
tures of the phase diagram, with the exception of the 
nematic-nematic coexistence, are similar to those pre- 
dicteiLby Lekkerkerker, Coulon, van der Haegen and De- 
bliektEI. The polydisperse Zwanzig model therefore pro- 
vides a useful starting point for understanding the effects 
of polydispersity on the phase behaviour of hard rod sys- 
tems. 

A further important motivation for studying polydis- 
persity, in this simplified model, is that it provides an in- 
teresting scenario for testing and extending the recently 
proposed "moment method" approach tO|-the-|thermody- 
namic treatment of polydisperse systemst3i23. The mo- 
ment method applies to systems whose excess free energy 
depends only on some moments of the density distribu- 
tion describing a polydisperse system; we show below 
that the Zwanzig model (treated within the second virial 
approximation) is of exactly this form. By expressing the 
ideal part of the free energy in a similar form a "moment 
free energy" can be defined, which only depends on the 
given moment densities. This is a drastic reduction in 
the number of densities required to describe the system, 
from the infinite number of degrees of freedom of the com- 
plete density distribution to a finite number of moment 
densities. The standard methods of the thermodynam- 
ics of finite mixtures can then be applied to the moment 
free energy to analyse the phase behaviour. Although in 
general the results will be approximate, for systems with 
conserved densities it has been shownE3~E3 that the cloud 
and shadow points (which specify the onset of phase co- 
existence in polydisperse systems), spinodals, and critical 
points are all found exactly. In the case of the Zwanzig (or 
Onsager) model, however, one has both conserved (rod 
lengths) and non-conserved (rod orientations) degrees of 
freedom. We show below that the moment method can 
be extended to this kind of scenario, and that most of 
the above exactness statements carry over. We also as- 
sess the accuracy of the moment method in the region 
where it is not exact and discuss how it can be improved 
systematically by retaining additional moment densities 
in the description. 

The paper is structured as follows. In Sec. ||, we de- 
scribe the polydisperse Zwanzig model, give its free en- 
ergy (calculated within the second virial approximation) 
and derive the corresponding expressions for the chem- 
ical potentials and the osmotic pressure. These results 



are used in Sec. Ill to study the exact thermodynamic 
behaviour of the model. In Sec. IV, we construct the 
moment free energy, discuss its properties, and compare 
the results obtained from it with the exact ones. Be- 
cause the exact calculation of the phase behaviour of the 
polydisperse Zwanzig model is feasible, it might seem 
unnecessary to confirm these results using the moment 
method. However, it is precisely because the exact re- 
sults are available that the Zwanzig model provides a 
useful test case for the application of the moment method 
to systems with non-conserved degrees of freedom. The 
conclusions drawn should help us to apply the method 
to more complicated systems (such as the polydisperse 
Onsager model) where an exact calculation of the phase 
behaviour is infeasible. This and other possible avenues 
for future work are discussed in Sec. 



II. DEFINITION OF THE MODEL 

The Zwanzig model considers hard rods in the shape 
of parallelepipeds of length L and with square base of 
edge length D. In contrast to the full Onsager model, 
the orientations of these rods are restricted to be along 
one of the three cartesian coordinate axes, x, y or z. We 
will assume that all rods have the same diameter I?, but 
that they are polydisperse in length, so that there is a 
continuous distribution of rod lengths L. Introducing a 
reference length Lq and the normalized lengths I = L/ Lq, 
we will focus on the Onsager limit of long thin rods. This 
corresponds to letting D/Lo while keeping the nor- 
malized lengths I constant. Unless a distinction between 
normalized (/) and unnormalized (L) lengths needs to be 
made explicitly, we will simply refer to / as the length of 
a rod in the following. 

Because of the length polydispersity, the number den- 
sities of the rods in the three possible orientations are 
specified by density distributions Pd{l), with ci = x, y, z; 
for a small range of lengths dl, Pd{l) dl is the number den- 
sity of rods oriented along d and with lengths between I 
and I + dl. The total number density distribution (irre- 
spective of orientation) is then 

p{i) = p^{i) + py{i) + p,{i) = j2 Pd{i) 



and integrating over I gives the total number density of 
rods: 



p= dl 



(1) 



Here and in the following, all integrals over I run from 
to oo. With these definitions, the density distribution 
Pd{l) over lengths / and orientations d can be decomposed 
as 



Pdil) = pil)Pi{d) = pP{l)Pi{d) 



(2) 
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where 



P{l)=p{l)lp 



is the normaUzed distribution of rod lengths I, and 



Pi{d) 



P{1) 



is the probabihty of finding a rod with given length I 
in orientation d. Note that Pi{d) is the analogue of the 
orientation distribution in the full Onsager model and 
obeys the normalization 
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(3) 



In the isotropic phase, Pi{d) = 1/3 for all d and I. In 
the nematic phase, on the other hand, we have Pi(x) — 
P;(y) < Pi(z) if we take the director to be along the 
z-axis. We will nevertheless develop the theory of the 
model first for arbitrary orientation distributions (with 
P;(x) ^ Pi{y)) because this leads to somewhat more com- 
pact expressions, and only specialize to the nematic case 
at a later stage. 

For simplicity, we only treat the model in the second 
virial approximation. (In contrast to the case of the On- 
sager model, this approximation does not become exact 
in the limit D^Lq — s- here: higher order virial terms 
to not vanishlij.) The excess free energy is then essen- 
tially determined by the excluded volume of two rods. 
If the rods have (normalized) lengths I and V and are 
perpendicular, this volume is 

Vl'^'^i = 2D{L + D){L' + D)^ 2ll'{DLl)[l + 0{D/Lo)] 
while the excluded volume for parallel rods 



y^cxcl 



D 



iD^L+L') = A{l+l'){DLl)— = DLlO{D/Lo) 



La 



0. As 



we choose DLq as our 



is negligible by comparison in the limit D/Lq 
suggested by the result for V^^'^^ 
unit of volume in the following, making all densities di- 
mensionless {p pDLq). If we also set ksT — 1, the 
excess free energy density becomes (within the second 
virial approximation) 



/ = 2 jdldl'll'[p^{l)Py{l')+P^{l)pS')+Py{l)pS')] 

= 2(0x0y + + 4>y4>2) 

(4) 

d 

where we have defined 

c^d^ jdllpdil)^ ldllp{l)Pi[d) (5a) 

(5b) 



Y,<t>d= jdllpil) 



Our choice of volume units implies that {D/Lo)4>d and 
{D/LQ)(j) are, respectively, the volume fraction of rods 
pointing in direction d and the total volume fraction. 
The ratio 



m = <j)/p = jdllP{l) 

is then simply the average length of rods in the system, 
i.e., the first moment of the rod length distribution P{1). 

Adding the ideal part of the free energy (density) to 
Eq. (^), we have for the total free energy (density) 



/ = E fdiPdm^Pdii)-i]+Y.' 

.1 -J .1 



M (6) 



This equation is the starting point of our analysis. When 
we refer to "exact" results in the following, we mean the 
exact thermodynamics of the model defined by the free 
energy ( 3|) . 

Eq. (p) shows that the free energy is a functional of the 
density distribution pd{l) over / and d. If we use Eq. (|^) to 
write pd{l) = pil)Piid), we note an important difference 
between the two factors: while the total density distri- 
bution p{l) is conserved (because the rods cannot change 
length) , the orientation distribution Pi (d) is not (because 
the rods can change orientation). We separate out the 
respective contributions to the free energy by writing 

/ = Jdl p{l) [In p{l) - 1] + jdl p{l) E Pi {d) In Pi (d) 

d 

+ Y,M^-<l>d) (7) 

d 

For a given p{l), the orientation distributions Pi{d) (for 
each I) are then obtained by minimizing /, subject to 
the normalization constraints (^) (again, for each I). In- 
troducing Lagrange multipliers k{1) for these constraints 
gives the minimization condition 



/ + 



lduii)Y,Piidyj = 



SPiid) 

p(Z)[lnP,(d) + 1] + 2(0 - cj)dMl) + k(0 = 

Solving for Pi{d) and eliminating the k{1) by using Eq. 
gives the orientation distributions 



Piid) = 



o2((/)d-0)/ 



(8) 



Inserting this result into the definition (Q) then gives 
three simultaneous nonlinear equations which can be 
solved for the (pd- 

To derive the conditions for phase coexistence in the 
polydisperse Zwanzig model, we need expressions for the 
chemical potential p{l) — which, due to the polydisper- 
sity, is a function of the rod length I — and the osmotic 
pressure. The chemical potential 



3 



m(0 = 



Sf 
5p{l) 



is obtained by functional differentiation of the free en- 
ergy (0) w.r.t. p{l). There is no contribution from the 
variation of Pi{d) with p{l) because we have aheady min- 
imized the free energy w.r.t. Pi{d). This leads to 



p{l) = \np{l) + ^'(^) ^^Piid) + Sj]!'^ - (^d)Pi{d)l 



or, using Eq. (|), 



= \np{l) - In e2('^''-^)'j 



(10) 



The osmotic pressure can be written in terms of the free 
energy and the chemical potential; hence, using Eq. (||), 



n = -/ + Jdip{i)p{i) = p + J2M 



(11) 



III. EXACT PHASE COEXISTENCE 
CALCULATION 

A. Coexistence conditions 

We can now state the conditions for coexistence of two 
or more phases, labelled by a = 1 . . . K, into which a 
"parent" phase with density distribution p''^\l) is as- 
sumed to have split. From Eq. (|ic|), the equality of the 
chemical potentials between the phases is obeyed exactly 
if the densities can be written in the form 



i?(0^exp 



(12) 



with a function R{1) common to all phases and the a^""* 
obeying 



)+c 



(13) 



Here c is an arbitrary constant (again common to all 
phases). If the phases occupy fractions v^°-^ of the total 
system volume, particle conservation implies 



Y^v^'^^p(^Hi)-P^"Hi) 



(14) 



The density distributions over rod lengths I and orienta- 
tions d are then found from Eq. §) and, using Eqs. (p|,[I3[), 
take the simple form 



p(°)(0 



exp 








d' exp 





(15) 



Integrals over these distributions define, by Eqs. (|]j|), 
the values of the densities p*^"^ and volume fractions (pj^^ , 

(9) 

(j)^"-' in all phases. These variables determine the pres- 
sures 



(16) 



in the different phases; at phase coexistence, these must 
of course all be equal. We thus have, in the most general 
form of the conditions for coexistence of K phases, AK 

variables (three a^^^ and one w*^"-' per phase a — 1 . . . K) 
and equally many equations to solve: the iK condi- 
tions (O) f or chemical potential equality, the K —\ con- 
ditions (p-6[) for equality of the pressures, and the trivial 
normalization of the phase volume fractions, v^°'^ = 1. 

It is easy to show that, just as in the Onsager 
model, isotropic-isotropic coexistence is not possible in 
the Zwanzig model with length polydispersity only. (This 
would be different if the rod diameters were polydisperse 
as well; compaxe Refs. pl|-p^.) Given the results for the 
bidisperse casetll, it is also unlikely that nematic-nematic 
coexistence could occur; this is in_contrast to what has 
been found for the Onsager modelEEl. Intuitively, the dif- 
ference can be explained as follows: when a polydisperse 
nematic phase splits into two nematics containing pre- 
dominantly short and long rods, respectively, it gives up 
entropy of mixing but gains orientational entropy. In 
the Onsager model, where the rod angles are continu- 
ous variables, the gain in orientational entropy can be 
arbitrarily large, thus favouring such a phase split. (The 
orientational entropy tends to — oo as the orientational 
distribution function tends to a delta function.) In the 
Zwanzig case, on the other hand, the maximum gain in 
orientational entropy is fee In 3 (this being the difference 
between the entropies of an isotropic and a fully ordered 
nematic phase) so that nematic-nematic coexistence is 
disfavoured. 

We therefore now specialize to coexistence between an 
isotropic (I) and a nematic (N) phase. If we choose the 
director to be along the z-axis, we then have (/)x = 4>y 
and = 2(/)x + </>z- Denoting 



This fixes R{1), giving 



p^-\l)^p^^'^{l) 



io). 



Ed exp 






Ea' "^"'^ Ed' exp 





the volume fractions of rods with the three possible ori- 
entations can be expressed as 



2A) (17) 
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and the excess free energy and pressure simplify to 



- A2) 



(18) 
(19) 



Instead of numbering the phases by a — 1,2, we label 
them with superscripts I and N from now on. In the 
isotropic phase, we have = <p^ /?> for d — x,y,z (A^ — 0), 
and by choosing the arbitary constant c in Eq. (0) as c = 
4(/)i/3 we can ensure that all the coefficients a\ vanish. If, 
for the corresponding coefficients in the ncmatic phase, 
we write 



the conditions 



N _ 
a = q;j_. 



3|) simplify to 



N _ 

a, = ail 



II 3^ 



+ A) 

--A) 
2 ' 



(20) 
(21) 



The condition ( p^ ) of equality of the pressures, on the 
other hand, becomes 



N^2 



A^ 



(22) 



Note that we have dropped the subscript 'N' on a||, 
a±_ and A because the corresponding quantities in the 
isotropic phase are all zero. 

If we denote = 7 and = 1 — 7 (so that 7 is the 
fraction of the system volume occupied by the isotropic 
phase), our phase coexistence problem now takes the 
form of three nonlinear equations ( |2C| , |2l| , ^ ) for a||, aj^ 
and 7. The densities and volume fractions appeari ng in 
these equations can be found by specializing Eqs. (p],p|]r^) 
to the case of I-N coexistence: 



with 



P\l) 



p'^ jdlp'il) 


(23a) 


(j)^ = Jdllp\l) 


(23b) 


= Jdlp^il) 


(23c) 


cb"" = Jdllp^'il) 


(23d) 


A = Idll[p]}{l)-p1il)] 


(23e) 


3pW(0 

^ 37 + (1-7) [e"ii' +2e°^'] 


(24a) 


= Jp^(Oe""' 


(24b) 


= ^P^(Oe"-' 


(24c) 




(24d) 



In the above setting of the phase coexistence problem, 
we specified a single parent density distribution /9^''^(Z). 
We will normally be interested in results along a so- 
called "dilution line", where the parent length distribu- 
tion F(°)(/) = (;)//?('') is kept fixed while the overall 
parent density p^^^ is varied. As p^^^ is increased from 
zero, we then expect to find a single isotropic phase first 
(7 = 1). At the isotropic "cloud point", an infinitesimal 
fraction of nematic phase will first appear; the density 
p^ of this nematic phase gives the nematic "shadow". 
On the other hand, starting from high density p^^^ we 
will first see a pure nematic phase (7 = 0). On de- 
creasing p*^'^) , an infinitesimal amount of isotropic phase 
will appear at the nematic cloud point; the density of 
the isotropic phase at this point gives the corresponding 
isotropic shadow. The two cloud points delimit the co- 
existence region. For values of p^^^ inside this region, an 
isotropic and a nematic phase coexist and occupy non- 
infinitesimal fractions of system volume, with 7 decreas- 
ing from 1 to as p^^-* increases. 

Numerically, rather than changing p(°) and finding 7, 
it is easier to vary 7 between and 1 and find the corre- 
sponding p^*^\ To implement this scheme, one only has 



(o)(;) = p(o)p(o)(;) 

and p'^°\ Alter- 
as being defined 



to replace p^ '{1) in Eqs. (gj) by p 
and solve Eqs. (|^,|l|,^ for a a 
natively, one can interpret a\\ and a 
by Eqs. (|2|,|2^), and regard p(°), p^ 0\ p^, 0^ and A 
as the underlying variables. Eqs. (^2|j2^) then constitute 
six equations for these six unknowns, which can be solved 
numerically; this is the approach that we adopt. 



B. Results: Cloud point and shadow curves 

In the following, we will restrict ourselves to the case 
where the rod lengths in the parent phase are distributed 
according to a Schulz distribution 



p(")(0 



2 + 1 



r(z + i) 



Pexp [~{z + l)l] 



(25) 



This distribution is normalized and has an average rod 
length of m^^^ — 1. (Allowing other values of m(°) would 
not make our treatment more general since the value of 
to'^*'^ can always be absorbed into a rescaling of the ref- 
erence length Lq.) The parameter z controls the shape 
and width of the distribution, and is taken to be non- 
negative. A more intuitive measure of the width of the 
parent distribution is the relative standard deviation a 
(usually called the "polydispersity"), defined by 



1 



7,(0)12 



1 



(26) 



It is then easy to see that, for the Schulz distribution, 
a=(l + z)-i/2 



5 





FIG. 1. The isotropic and nematic cloud point curves 
(solid) and their corresponding shadow curves (dashed), show- 
ing the densities of the coexisting phases as a function of the 
polydispersity a. Of the two cloud point curves, the isotropic 
one is the one with the lower density; it meets the isotropic 
shadow curve for ct as it must. The nematic cloud and 
shadow curves likewise coincide in this limit. 



FIG. 2. The average rod lengths m in the shadow phases 
are plotted as a function of the polydispersity, o. The dashed 
and solid curves give the results for the nematic and isotropic 
shadow phases, respectively. Note that the corresonding cloud 
phases are identical to the parent and therefore have average 
rod lengths equal to m'"-* = 1. 



For z ^ CX3, wc thus have a monodisperse parent with 
CT = and p(")(/) = 5{l — 1). As z is decreased, a in- 
creases and the parent gets more and more polydisperse. 
For = 0, finally, the parent distribution is a simple 
exponential, and a achieves its (for the chosen Schulz 
distribution) maximal value of 1. 

To calculate the cloud poin t and shadow curves, we 
proceed as explained in Sec. 



Ill A 



For the isotropic 
cloud point and shadow, we set 7 ~ 1. In Eqs. (24), 
we then have p^{l) — p''^\l). This of course makes sense: 
only an infinitesimal amount of nematic phase has ap- 
peared, and so the density distribution of the isotropic 
phase is only negligibly perturbed away from the parent. 
In Eqs. (p3|), the equations for and (f>^ then simplify to 
the trivial statements p^ — p^"^ and 0^ — p'^^^nS'^^ — p'^^\ 
and we only have to solve four equations for the four 
unknowns p^'^\ p^ , (j)^ and A'^. Conversely, for the ne- 
matic cloud point and shadow, we set 7 = 0. We then 
find p^{l) = p^^^l) and p^ = (j)^ = p^^^ and have to solve 
the remaining four equations for the four unknowns p^^^ , 
p\ 0i and A^. 

The results for the cloud point and shadow curves are 
plotted in Figs. |l| to ^. In Fig. |^ we see that the co- 
existence region (the density range between the isotropic 
and nematic cloud points) broadens quite dramatically as 
the parent distribution becomes more polydisperse. The 
transition, which is already strongly first order in the 
monodisperse case, spreads out so that when a = 1 the 
coexistence region spans almost an order of magnitude in 
density, from p — 0.54 to p = 3.96. As a increases, the 
nematic shadow curve moves rapidly towards lower densi- 
ties, approaching the isotropic cloud curve. In Fig. ||, we 



show the average rod lengths in the nematic and isotropic 
shadow phases. As the polydispersity increases, a strong 
fractionation effect is observed, with long rods found 
preferentially in the nematic phase; this is in qualitative 
agreemeat with results for the bi- and tridisperse Onsager 
modeinl3. For cr = 1, for example, the average length of 
rods in the nematic phase, m^, is more than double that 
in the isotropic phase, m}, both at the isotropic and at 
the nematic cloud point. (Note that at the isotropic cloud 
point, m} = TO^"-* = 1 and > 1, while at the nematic 
cloud point, = m^^^ = 1 and < 1-) This fraction- 
ation effect can be seen in more detail in Fig. ||, where 
for a — 0.75 we have plotted the relevant rod length 
distributions P{1). For example, a,t I — 4 {i.e., at four 
times the average rod length of the parent), P{1) in the 
nematic shadow phase is almost an order of magnitude 
larger than in the isotropic cloud phase. Finally, we 
study in Fig. ^ the cloud point and shadow curves in a 
different representationcJ: instead of the number density 
p of the coexisting phases (as in Fig. |l|), we show their 
rescaled rod volume fraction (j) = mp. This leaves the 
cloud point curves (for which m = m^"^ = 1) unchanged, 
but does affect the shadow curves (along which, as Fig. || 
shows, TO can differ significantly from m^^^ = 1). Inter- 
estingly, the volume fractions of the shadow phases turn 
out to depend only weakly on the polydispersity a, in 
contrast to their number densities (compare Fig. In 
fact, the volume fraction in the nematic shadow phase 
does not even show a definite trend in its dependence on 
polydispersity, being a non-monotonic function of a. 

We have not so far discussed the strength of the ori- 
entational ordering in the nematic phase. This can be 
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FIG. 3. The rod length distributions P{1) in the coexisting 
isotropic and nematic phases, at polydispersity a = 0.75 and 
for different fractions 7 of system volume occupied by the 
isotropic phase. Bold dashed line: P{1) in the nematic shadow 
at the isotropic cloud point, 7 = 1. The distribution in the 
isotropic phase at this point is the parental one (bold solid 
line). As 7 decreases through 0.75,0.5,0.25, the dashed and 
dotted lines show — from bottom to top — the distributions in 
the coexisting nematic and isotropic phases, respectively. At 
7 = 0, finally (the nematic cloud point), the nematic has the 
parental length distribution; the bold dotted line shows P{1) 
in the corresponding isotropic shadow. Inset: Ratio of the 
rod length distributions to that of the parent. 




FIG. 4. The isotropic and nematic cloud point curves 
(solid) and their corresponding shadow curves (dashed) of 



id) 

Fig. nl are replotted here, showing the (rescaled) rod volume 
fractions (p = mp of the coexisting phases rather than their 
densities p. Of the two cloud point curves, the isotropic one is 
the one with the lower (p; it meets the isotropic shadow curve 
for (T — > as it must. The nematic cloud and shadow curves 
likewise coincide in this limit. 



characterized by the order parameter q = A/(f>, which 
has the value g = in the isotropic phase and 5=1 for 
perfect nematic ordering. It turns out that, due to the 
strongly first order nature of the I-N transition, q is close 
to its maximal value of unity for all polydispersities cr, so 
we do not display it. 



C. Results: Inside the coexistence region 

We now turn to the properties of the isotropic and 
nematic phases in the coexistence region, i.e., for parent 
densities between the isotropic and nematic cloud points. 
Both phases then exist in non-infinitcsimal amounts, im- 
plying < 7 < 1. Using the numerical scheme out- 



lined in Sec. Ill A, we then obtain the results shown in 
Figs. H and ^. Fig. || tracks the densities of the coexisting 
phases as the coexistence region is crossed (for a = 0.5). 
As expected, the densities interpolate between the cloud 
and shadow phase densities at either end and increase 
smoothly with the parent density. Fig. |^ shows similarly 
the variation of the rod lengths in the isotropic and ne- 
matic phases across the coexistence region. As expected 
from Fig. H, the average rod length in the nematic phase, 
TO^, is always higher than that in the isotropic phase, 
m^. At the isotropic cloud point, = 1 and > 1, 
while at the nematic cloud point, m} < 1 and = 1; 
again the values inside the coexistence region interpolate 
smoothly between these limits, with both average rod 
lengths decreasing as the parent density p^"-* increases. 

Finally, we can also study the evolution of the distri- 
bution functions P{1) as the fraction of volume occupied 
by the isotropic phase, 7, is varied. The results are in- 
cluded in Fig. ||. At the isotropic cloud point (7 = 1) the 
isotropic phase has the parental distribution P*^°)(/); as 
7 is decreased (corresponding to increasing parent den- 
sity p'-^-'), this distribution shifts towards smaller lengths, 
evolving smoothly into the distribution at the nematic 
cloud point 7 = 0. Proceeding in the reverse direction, 
the nematic phase has the parent distribution at 7 = 
and then changes smoothly into the distribution at the 
isotropic cloud point as 7 is increased towards 1, shifting 
towards larger rod lengths in the process. 



IV. COMPARISON WITH THE MOMENT 
METHOD 

A. Constructing the moment free energy 

We now outline how the moment methodEi~0 can be 
applied to the polydisperse Zwanzig model. To construct 
the moment free energy, one recognizes from Eq. (|l^ ) 
that the excess free energy of the model (specialized to 
isotropic or nematic orientational order) only depends on 
the variables 6 and A. Both of these are moments of the 
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FIG. 5. The densities p of the coexisting isotropic (solid) 
and nematic (dashed) phases as a function of the parent den- 



sity p 



(0) 



for polydispersity a = 0.5. The isotropic and ne- 



matic cloud points, which delimit the coexistence region, are 
located at the densities where and p^ meet the "dilution 
line" p = p^°^ (dotted), respectively. Outside the coexistence 
region, there is only a single isotropic (for low densities) or 
nematic (for high densities) phase with density distribution 
p^'^\l) (and therefore density p'"') identical to that of the 
parent. 
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FIG. 6. The average rod lengths in the coexisting isotropic 
{m}, solid) and nematic {m^ , dashed) phases corresponding 
to Fig. p. The dotted line indicates the average rod length 
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= 1 of the parent. 



4>=Pi='^ jdlwi{l,d)pd{l) 
d •' 

A = p2 = 51 dlw2{l,d)pd{l) 

d 

defined by the weight functions 
wi(/, d) ~ I 

W2{1, d) = I (Sd,2. - -^Sd,x - 2^"*'^ 



We therefore call pi and p2 moment densities. Another, 
trivial, example of a moment density is the total number 
density 



which corresponds to the weight function wo{l,d) — 1. 
Even though the excess free energy 

depends on moment densities only, the ideal part 

/ideal = ^ dl pd{l)[ln pd{l) ~ 1] 
d •' 

of the free energy / = /ideal + / still contains all details 
of the density distribution Pd{l)- To construct a moment 
free energy which depends only on the moment densities 
appearing in /, we therefore need to transform this ideal 
part to a moment form. For this purpose, it is useful 
to add a term — Jdl pd{l) In r(/) — — Jdl p{l) In r(^) to 
the free energy, giving 



r{l) 



f 



The additional term is linear in the conserved densities 
p{l) and therefore has no effect on the exact thermody- 
namics described by /. (This would not be true if we 
had replaced r{l) by a d-dependent quantity rd{l), be- 
cause pd{l) is not conserved.) For the moment method, 
on the other hand, r{l) turns out to be crucial. The 
key idea is to allow violations of the particle conserva- 
tion rule ( p^ ) as long as they do not affect the moment 
densities appearing in the excess free energy (0 = pi and 
A = p2, in our case). The intuitive rationale — to be 
verified a posteriori — is that phase behaviour is mainly 
governed by the excess free energy and the moment den- 
sities appearing in it. The relevant free energy is then 
obtained by minimizing / at given values of the moment 
densities pi {i = 1,2). Using Lagrange multipliers \i to 
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fix these values, one finds that the minimum of / occurs 
for density distributions of the form 



= r{l) exp 



^\iWi{l,d) 



(27) 



and the corresponding minimum value is 

/mom = ^ \Pi - Po + / (28) 



This expression defines the moment free energy. Note 
that, from Eq. (^), the moment densities are related to 
the Lagrange multipliers by 



Pt = 



= dlwi{l,d)r{l)exTp 



Inverting these relations determines the Xj in terms of 



the Pi; the moment free energy ( ^8]) thus depends on the 
moment densities only, as desired. 

By construction, the moment free energy (^ ) is the 
free energy of systems with density distributions of the 
form (27). One therefore expects it to give exact results 
for the phase behaviour as long as the density distribu- 
tions of all the coexisting phases are actually contained 
in the "family" (p7|). Considering an isotropic parent, 

with pI^\i) = p'*'^(l)/3, we can ensure that this is true 
at least for the parent by choosing r{l) = p^^\l)/3. This 
is the choice that we adopt from now on, giving explicitly 



XiWi{l, d) 

i 

for the family of density distributions and 



(29) 



Jdlw,{l,d)^p^°\l)exp 



^XjWjild) 



(30) 



for the relation between the Lagrange multipliers Xj and 
the moment densities pi. With this choice, the isotropic 
cloud point and corresponding nematic shadow will be 
found exactly by the moment method: at that point, 
the parent is only negligibly perturbed because only an 
infinitesimal amount of the nematic phase has appeared, 
while the nematic phase itself is related to the parent by 
exactly the kind of Gibbs-Boltzmann factor appearing in 
Eq. (§). 

To see more formally why the moment method gives 
exact results for the isotropic cloud point, let us write 
down the resultant phase coexistence conditions and 
show that they are equivalent to the exact conditions 
(particle conservation violations apart). Associated with 
each of the moment densities pi is a moment chemical po- 
tential iJ,i p=j dfrnom/dpi- Using the Legendre transform 
propertiesE3 of /mom, one finds 



9f ^ 4 
Ml = Ai + - — = Ai + -pi 
dpi 3 

df X 4 
X2 + t; — = A2 - -P2 
op2 3 



(31a) 
(31b) 



Because pi is conserved while p2 is not, pi must be equal 
in all coexisting phases, while p2 must actually be zero: 



(a) (a) „ 



for all a 



(32) 



where c is a constant common to all phases. To com- 
pare these conditions with the exact conditions (12|T^) 
for equality of the chemical potentials p{l), we write the 
density distributions ( p9| ) for the different rod orienta- 
tions explicitly: 

p„(O^Pz(0 = ^p(°HOe(^^+'^" 

p^il)^p^{l)=p,il)^^p('\l)e('^~'^/'^' 
These are of the form (|lj) if we identify 



a I = ttx 



Ai - A2/2, ail = Qfz = Ai + A2 (33) 



Using Eq. ( p^ ) , coexisting phases calculated from the mo- 
ment free energy therefore have equal (exact) chemical 
potentials p{l) if 

a1°)+A^°) =2 ((/.('') ~0('^)) + c 

^ + AW) +c 

3 

A('^)-iA(")=2(4'^)-0("))H-c 

= -(-2(/i('^)-A('^))+c 
3 

in all phases. But from Eqs. (|3^) one easily sees that 
these conditions are equivalent to those [Eq. (|3^)] de- 
rived from the moment free energy, as promised. For 
the condition of equality of pressure in all phases, it is 
even easier to see that the moment free energy gives the 
correct answer: one finds 



n = -/„ 



PiPi + P2P2 = Po 



^{pI-pD 



which is exactly the same as the result (19) derived from 
the original free energy /. Remember that, in our mo- 
ment density notation, po = P, Pi = (t>, and p2 = A. 

We have thus shown that coexisting phases calculated 
from the moment free energy satisfy the exact phase equi- 
librium conditions of equal chemical potentials and pres- 
sures. If the phases also obey the exact particle conser- 
vation conditions (0), they therefore give the exact so- 
lution of the phase coexistence problem. This is the case 
at the isotropic cloud point, because one of the phases is 
then identical to the isotropic parent, and the other (ne- 
matic) phase is infinitesimally small. As stated above. 



9 



this point will therefore be located exactly by the mo- 
ment free energy method. The nematic cloud point, on 
the other hand, will not be found exactly: on the high- 
density side of this point, the density distribution of the 
single nematic phase — p^^\l)Pi{d), with Pi{d) 

obeying Eq. (||)), will not in general be a member of the 
family dH)- 

In Refs. |l8|,|9| it was shown that the moment free en- 
ergy allows one to determine exactly the onset of phase 
coexistence (cloud point and shadow), the spinodals and 
the critical points of a polydisperse system with con- 
served densities. In our above discussion, we have shown 
that for a system with non-conserved degrees of freedom 
(the rod orientations), the onset of phase coexistence is 
still located exactly under the following condition: in the 
single phase region from which coexistence is approached, 
the parent must not exhibit any ordering of the non- 
conserved degrees of freedoms (which means in our case 
that it must be isotropic). One can show that this con- 
clusion holds quite generally, and that under the same 
restriction spinodals and critical points found from the 
moment free energy also remain exact. We thus con- 
clude that the moment method remains useful even in 
systems with non-conserved degrees of freedom. 

To improve the accuracy of the moment method in the 
regions where it is not exact (beyond the isotropic cloud 
point), one can simply retain additional moment densi- 
ties in the moment description, defined by weight func- 
tions Wi{l,d). The above construction of the moment 
free energy generalizes directly to this case: the expres- 
sions ( p8|p9|j30| ) remain valid as long as the sums over 
i are extended appropriately. From Eq. (|29|), one rec- 
ognizes that the addition of new moment densities has 
the effect of extending the family of density distributions 
that are accessible; the exact distributions [Eq. (|l^)] in 
the coexisting phases can thus be approximated with ar- 
bitrary accuracy as the number of moment densities is 
increased. More explicitly, this can be seen as follows. 
Because the additional moment densities pi (i 7^ 1, 2) do 
not appear in the excess free energy, their associated mo- 
ment chemical potentials are simply pi — Xi. These must 
be equal in all phases, so we can write the density dis- 
tributions in coexisting phases predicted by the moment 
method as 



XiWi{l,d) 



X exp 



Comparing with Eq. (|12|), and bearing in mind the identi- 
fication (p3|), we see that the moment method essentially 
approximates ln(3/£)) (where D is the denominator on 
the r.h.s. of Eq. (|l^)) by a linear combination of the ad- 
ditional weight functions. Since D depends on / only 
(not on d), all these weight functions can be chosen to be 



d-independent; the corresponding moment densities are 
thus conserved. A particular additional moment density 
that we will always retain is the overall density po, with 
weight function wq{1) — 1. This guarantees that the dilu- 
tion line p{l) — const x p^^'^H) — e'^°p^^''{l) for the parent 
p("^(Z) is contained in the family (|27|), and thus simplifies 
the calculation of cloud points and shadows. The optimal 
choice of the remaining additional weight functions Wi{l) 
{i > 3) is less clear cut and is discussed further below. 
One thing we can say already at this point, however, con- 
cerns the large I asymptotics: it is easy to see that, for 
large I, ln(3/£>) = ci + C2I + e~'^^\ with constants ci, C2, 
C3, and up to terms which are exponentially smaller. The 
first two contributions are covered by by the weight func- 
tions wq, wi, W2, so all other weight functions should be 
chosen to decay exponentially for large I; the coefficient 
C3 of this decay is not known a priori, however. 



B. Results 

We show in Fig. |^ the cloud point and shadow curves 
obtained from the moment free energy with different 
numbers n of moment densities retained. As explained 
above we always keep, beyond the "essential" moment 
densities pi and p2, the overall number density po, so the 
smallest value of n that we consider is n = 3. For larger n, 
the additional weight functions were chosen to be expo- 
nentials with increasing decay constants, W2+j{l) — e~'^^^ 
(j = 1 ... n — 3). This form is consistent with the ex- 
pected exponential behaviour for large Z; the coefficient 
c was chosen as c = 0.5 by trial and error. As expected, 
the isotropic cloud point curve and the corresponding ne- 
matic shadow are found exactly already for = 3. For 
the nematic cloud point curve and the isotropic shadow, 
deviations from the exact results become apparent for 
larger polydispersities cr; as expected, these deviations 
decrease as n, the number of moment densities retained, 
increases. 

As noted in the introduction, for the polydisperse 
Zwanzig model that we are considering the moment 
method is not really needed because the exact phase co- 
existence equations can be solved directly. However, for 
more realistic models (such as the polydisperse Onsager 
model, with unrestricted rod orientations), this will not 
be the case. One thus needs to be able to assess the 
accuracy of the moment method without knowing the ex- 
act results beforehand. In Ref. the following quan- 
tity was proposed for this purpose: for any phase coex- 
istence calculated from the moment free energy, one can 
work out the total density distribution over rod lengths, 

Ptot(0 = Ea^^^V^-^HO = Ea,d^^'^Vi"^(0- The quan- 
tity lnptot(0/P^°''(0 then measures how strongly particle 
conservation of rods of length / is violated; taking the 
square and averaging over the normalized parent distri- 
bution P(o)(0 = p(o)(0/p(°) defines the "log-error" 
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FIG. 7. Cloud point and shadow curves found using the 
moment method. The exact results (compare Fig. ^) are 
shown in bold for comparison. Even when only the mini- 
mal number of moment densities (n = 3) is retained in the 
moment free energy, the isotropic cloud point and correspond- 
ing nematic shadow are found exactly; the moment method 
results therefore overlay the corresponding exact curves. The 
nematic cloud point and isotropic shadow are not found ex- 
actly, but their accuracy increases as n (indicated by the up- 
per row of symbols) increases. The other symbols show the 
points on the moment method curves where the log-error 5 
first reaches the value 10"* as a is increased from zero; see 
text for discussion. 



8 = jdlP^°\V) [\n 



Ptot(0 
p(0)(/) 



(34) 



For small violations of particle conservation, 
lnptot(0/p'"HO ~ /Otot(0/p'"'(0-l, and we can think of 
\/5 as the root-mean-squared relative deviation between 
ptot(0 and p^'^^l). In Fig. 0, we indicate by the lower 
symbols on the nematic cloud point and isotropic shadow 
curves where, as the polydispersity a is increased from 
zero, 6 first reaches the value 10~^ (which corresponds 
to an average violation of particle conservation of 1%). 
The fact that the symbols lie essentially on the curves 
with the exact result shows that 5 provides argpod indi- 
cator of the accuracy of the moment methodE3: adding 
moment densities until 6 < 10~^ ensures that the results 
are essentially indistinguishable from the exact ones. 

Beyond the calculation of phase boundaries, one would 
also like the moment method to give reliable results for 
the properties of coexisting phases inside the coexistence 
region. In Fig. ||, we therefore show the analogue of Fig. ^ 
for the moment method: the densities of the coexisting 
isotropic and nematic phases in the coexistence region, 
as a function of the parent density p^^\ Again, exact 
results are obtained only at the isotropic cloud point; but 
as the number of moment densities, n, is increased, the 



FIG. 8. The densities p of the coexisting isotropic (solid) 
and nematic (dashed) phases as a function of the parent den- 
sity p'"', for polydispersity a = 0.5. We show here the re- 
sults calculated from the moment method with n = 3,4,5 
(thin lines), and the exact results of Fig. ^ (bold lines). As 
expected, the densities beyond the isotropic cloud point are 
not exact, but become increasingly more accurate as n is in- 
creased. The results for n = 5 are indistinguishable from the 
exact ones on the scale of the plot. 



results across the whole of the coexistence region become 
progressively more accurate. 

We compare different choices for the additional weight 
functions in Fig. ^ in terms of the dependence of the 
log-error 5 on n at a point deep within the coexistence 
region of the phase diagram. Results for two sets of 
additional weight functions are shown. The first set 
consists of the exponential weight functions considered 
above. For the second set, we chose increasing powers of 
I, W2+j{l) = P-^e-"^ {j = l...n- 3). Note the expo- 
nential factor, which ensures the required asymptotic be- 
haviour; c was again chosen as 0.5. We call these weight 
functions "power-exponential" . 

In this example, the exponential weight functions are 
seen to lead to a faster decrease of 5 with n. However, fur- 
ther experimentation with other choices of weight func- 
tions may well lead to even better results, and it would 
clearly be desirable to have a more systematic way of 
constructing optimal additional weight functions. The 
following adaptive approach is a first step in this direc- 
(see also Ref. 



tion (see also Ref. 20). Consider a given point in the 
phase diagram, characterized in the present case by the 
density p'"^ of the parent and its polydispersity a. Per- 
forming a moment method phase equilibrium calculation 
without additional weight functions, one will find a cer- 
tain log- lever rule violation In ptot (0/p(°HO (called "log 
ratio" for short in the following) . One then expects that 
adding a weight function (w^{l), in our case) which has 
the same /-dependence as this log-ratio should signifi- 
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FIG. 9. The dependence of the log-error S on the number 
n of moment densities retained in the moment free energy, at 
the point p'"' = 2, cr = 1 in the phase diagram. The two lines 
correspond to different choices of the additional weight func- 
tions: exponential (solid) and power-exponential (dashed); 
see text for precise functional forms. 

cantly reduce the log-error: it extends the family ( p9| ) of 
density distributions "in the right direction" . Of course, 
constructing W3{1) to fit the log-ratio exactly would be 
computationally costly. Instead, we represent it as a lin- 
ear combination '^Ckipkil) of some simple basis func- 
tions ipk{l)- These could be the exponential or power- 
exponential functions used above, for example. The coef- 
ficients Cfc are chosen to minimize the squared deviation 
from the log-ratio (weighted by the normalized parent 
distribution) , 

This is a straightforward weighted least squares prob- 
lem, and the Ck can easily be found in closed form, thus 
determining w^^l). One can now repeat the phase equi- 
librium calculation with the moment defined by ^3(0 
included, and fit a new weight function 104(1) to the re- 
sulting log-ratio (which is expected to be rather smaller 
in magnitude than before). Repeating this process should 
lead to a steady decrease of the log-error 6. However, a 
large number of additional weight functions may still be 
required before S reaches an acceutably small value, and 
this can cause numerical problem£j. To avoid this prob- 
lem, we note from the discussion at the end of the previ- 
ous section that a single additional weight function can 
reproduce the exact results within the moment method, if 
only its l-dependence can be found appropriately. Rather 
than keeping a large number of additional weight func- 
tions, we can thus continually adapt a single weight func- 
tion, as follows. We choose the first additional weight 



function w^^l) by fitting the initial log-ratio, rerun the 
phase equilibrium calculation with included, and fit 
a "temporary" additional weight function W4{1) to the 
resulting decreased log-ratio. With both p^ and p^ in- 
cluded, we again run the calculation; this produces values 
of the Lagrange multipliers A3 and A4 (which, being asso- 
ciated with moment densities not appearing in the excess 
free energy, are common to all phases). The key point is 
now that if we merge w^ll) and w^^l) into the linear com- 
bination W3(/) = A3U'3(Z) -I- \4W4(l), and discard w^^l), 
repeating the calculation would give exactly the same re- 
sults. (All moment phase equilibrium conditions are still 
satisfied, and the lever rule for p'^ obviously follows from 
that for p3 and p^.) We are now back to a situation with 
only a single additional weight function and can repeat 
the process: obtain ^^(Z) by fitting to the current log- 
ratio, rerun with and p'^, combine w'^{l) and w'^ll) into 
w'^{l) and so on. This method avoids the computational 
problems associated with using a large number of addi- 
tional moment densities; indeed, it requires at most two 
additional weight functions at any time. The number 
of basis functions, however, is unrestricted in principle 
and can feasibly be made quite large. In fact, one can 
show that for an infinitely large set of basis functions, 
which allows arbitrary functional forms of the log-ratio 
to be fitted, the method must converge to the results of 
the exact phase equilibrium calculation (assuming it con- 
verges at all) . For finite but sufficiently large sets of basis 
functions, one thus expects excellent approximations to 
the exact results. As long as the set of basis functions 
is sufficiently "flexible" to approximate the /-dependence 
of the log-ratio, the precise choice of the basis functions 
should also be relatively unimportant, thus reducing the 
effect of the remaining heuristic element of the method. 

In Fig. we show the results for the adaptive method 
just described, at the same point in the phase diagram 
as in Fig. ||. As basis functions we considered the expo- 
nential and power-exponential weight functions described 
above. The number of basis functions was chosen such 
that if (in the previous non-adaptive, "brute-force" ap- 
proach) all basis functions are retained as additional 
weight functions, the log-error S is less than 10^^. As 
can be read off from Fig. Q, this leads to n — 3 = 7 — 3 = 4 
exponential basis functions and n — 3 = 10 — 3 = 7 
power-exponential basis functions. The corresponding 
"brute-force" values of S are shown as horizontal lines 
in Fig. |l^. They provide natural baselines for the results 
of the adaptive method: because the latter only retains a 
single additional weight function (a linear combination of 
the basis functions), it can obviously do no better than 
the brute-force method which allows the coefficients of 
all basis functions to be adjusted individually. Fig. |l^ 
confirms this; the adaptive method converges after a few 
iterations to a value of 5 above the brute-force baseline. 
While the slightly larger final value of 5 is, of course, 
a disadvantage, the adaptive method more than makes 
up for this by being much faster and numerically more 
stable. We therefore plan to study this method in more 
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FIG. 10. The dependence of the log-error S on the number 
of iterations of the "adaptive" weight function algorithm de- 
scribed in the text, at the point p'"' = 2, a = 1 in the phase 
diagram. Left: using four exponential basis functions; right: 
using seven power-exponential basis functions. Dotted lines: 
Value of 6 reached by the "brute-force" approach where all 
basis functions are retained as weight functions. 



detail in future. In particular, if one is interested in per- 
forming calculations for a number of points in the phase 
diagram falong a dilution line, for example, where the 
density p^^^ of the parent is varied) one could imagine 
a dynamical version of the algorithm which adapts the 
single additional weight function whenever a threshold 
value of the log-error S is crossed. As long as the chosen 
set of basis functions is sufficiently powerful, this should 
lead to uniformly precise results across the whole phase 
diagram. 

Finally, we illustrate in Fig. 
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the geometrical intu- 
ition provided by the moment free energy. We plot /mom 
as a function of the "essential" moment densities pi, p2, 
for a parent p^°\l) whose density was chosen to be ex- 
actly at the isotropic cloud point. As expected, the tan- 
gent plane drawn at the parent touches the surface at 
a second point: the nematic shadow phase. The mo- 
ment free energy thus allows a simple geometrical inter- 
pretation of this phase transition in a polydisperse sys- 
tem, in terms of a double-tangent plane to a conventional 
two-dimensional free energy surface. We emphasize that 
the properties of the cloud and shadow phases are found 
exactly, even though the moment free energy is only a 
low-dimensional projection of the true free energy (which 
"lives" in the infinite-dimensional space of density distri- 
butions Pd(l))- 



V. CONCLUSION 

We have studied the phase behaviour of the Zwanzig 
model of suspensions of hard rods, allowing for polydis- 
persity in the lengths of the rods. The model assumes 




FIG. 11. Moment free energy /mom versus the "essential" 
moments densities pi and p2, for a Schulz parent with a = 0.5 
and density p'"' = 0.959 corresponding to the isotropic cloud 
point. Constant and linear terms have been added to /mom to 
make its tangent plane at the parent (represented by the point 
pi = p'^'m'*^' — 0.959, p2 = 0) coincide with the xy-pl&ne. As 
expected, this tangent plane touches the free energy surface 
at a second point: the nematic shadow, whose values of pi 
and p2 are found exactly. 



that the rods are restricted to lie along one of three 
orthogonal axes. In spite of this drastic simplification 
(compared to the Onsager model, where rod orientations 



are unrestricted), the results we obtain are-in qualitative 
agreement with experimental observationsQEl: the coexis- 
tence region broadens significantly as the polydispersity 
(the width of the rod length distribution of the parent) 
increases; fractionation is also observed, with long rods 
found preferentially in the nematic phase. These conclu- 
sions were obtained from an exact analysis of the phase 
equilibrium equations, starting from the free energy of 
the model within the second virial approximation. 

In the second part of the paper, we considered the 
application of the moment method to the polydisperse 
Zwanzig model. This involved extending the construction 
of the moment free energy to models with both conserved 
and non-conserved degrees of freedom. We showed that 
most of the exactness statements obtained previously for 
systems with conserved densities carry over to this case: 
the onset of phase coexistence is still found exactly from 
the moment free energy, as long as it is approached from a 
single phase region where there is no ordering of the non- 
conserved degrees of freedom. With the same restriction, 
spinodal instabilities and critical points are also located 
exactly. Our concrete results for the cloud and shadow 
curves bear this out; the isotropic cloud point and cor- 
responding nematic shadow are found exactly, while the 
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nematic cloud point and isotropic shadow are approxi- 
mate. The accuracy of the approximation increases as 
the number of moment densities retained in the moment 
free energy is increased. The log-error (5 is a useful crite- 
rion for monitoring the increase in accuracy; crucially, it 
can be computed without knowing the exact results be- 
forehand. Finally, we have discussed methods for choos- 
ing the weight functions of the additional moment den- 
sities. An adaptive technique, which requires at most 
two additional weight functions at any given time, and is 
therefore very cheap to implement computationally, gives 
promising results. 

In future work, it is obviously desirable to remove the 
simplications of the Zwanzig model and move towards 
a study of the polydisperse Onsager model, with unre- 
stricted rod orientations. A direct numerical solution 
of the phase equilibrium equations for this model is in- 
feasible, so an approach based on the moment method 
suggests itself. One complication is that the excess free 
energy does not have a simple moment structure: the 
excluded volume between two rods at angles 9 and 0' 
with the nematic axis is a non-trivial function of thes€| 
angles. Its expansion in terms of Legendre polynomialsESI 
shows that the moment free energy actually depends on 
an infinite number of moment densities. As an interme- 
diate step, we thereiare plan to consider the polydisperse 
Maier-Saupe modelEZi, which truncates the eigenfunction 
expansion after the leading term and leads to an excess 
free energy depending on only two moment densities. 
This approach should be of independent interest as a phe- 
nomenological description of polydisperse suspensions of 
rod-like particles with more complex (soft) interactions. 
Finally, it would also be interesting to study the effect of 
diameter-polydispersity on hard rod systems. Novel fea- 
tures such as isotropic-isotropic phase coexistence have 
previously been found fQc-itiic bidisperse case (rods with 
two different diameters)EJ"L3, and it will be interesting 
to see how these are modified for truly polydisperse sys- 
tems. 
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